rm(list=ls())

set.seed(1234)

library(plyr)
library(ggplot2)
library(car)
library(psych)
library(texreg)
library(effectsize)
library(xtable)

setwd(".../Data/Analysis")

#read in country level datasets 

MexAgg<-read.csv("MexAgg.csv")
PeruAgg<-read.csv("PeruAgg.csv")
UKAgg<-read.csv("UKAgg.csv")
SpainAgg<-read.csv("SpainAgg.csv")
PortAgg<-read.csv("PortAgg.csv")
USAgg<-read.csv("USAgg.csv")
BrazilAgg<-read.csv("BrazilAgg.csv")
ArgAgg<-read.csv("ArgAgg.csv")
NZAgg<-read.csv("NZAgg.csv")
NorAusFranceAgg<-read.csv("NorAusFranceAgg.csv")

#make sure that names are the same across country level datasets to prepare for merging
names(UKAgg) <- names(MexAgg)
names(SpainAgg) <- names(MexAgg)
names(PortAgg) <- names(MexAgg)
names(BrazilAgg) <- names(MexAgg)
names(ArgAgg) <- names(MexAgg)
names(NZAgg) <- names(MexAgg)
names(NorAusFranceAgg) <- names(MexAgg)

#merge (rbind) all that sets 
MergedN<-rbind(MexAgg, PeruAgg, UKAgg, SpainAgg, PortAgg, USAgg, BrazilAgg, ArgAgg, NZAgg, NorAusFranceAgg) 

#ensure treatment numbers are hte same across datasets 
MergedN$treat<-recode(MergedN$treat, "  'treat1' = 1; 'treat2' = 2; 'treat3' = 3; 'treat4' = 4; 'treat5' = 5; 'treat6' =6")

#recode variables for consistency across datasets 

MergedN[MergedN =="Female"]<-1 #female coded as 1
MergedN[MergedN =="Feminino"]<-1 #female coded as 1
MergedN[MergedN =="Mujer"]<-1 #female coded as 1
MergedN[MergedN =="Woman"]<-1 #female coded as 1
MergedN[MergedN =="Hombre"]<-0 #female coded as 1
MergedN[MergedN =="Male"]<-0 #female coded as 1
MergedN[MergedN =="Man"]<-0 #female coded as 1
MergedN[MergedN =="Masculino"]<-0 #female coded as 1
MergedN[MergedN =="Other"]<-NA #female coded as 1
MergedN[MergedN =="Otro/ No binario"]<-NA #female coded as 1

#recode values on political ideology, for consistency across datasets 
MergedN[MergedN =="Left\n(0)"]<-0 #
MergedN[MergedN =="Right\n(10)"]<-10 

MergedN$leftright_1<-as.numeric(as.character(MergedN$leftright_1))

MergedN$gender<-as.numeric(as.character(MergedN$gender))

#subset dataset for those who saw a vignette about sexual harassment 
harass<-subset(MergedN, treat == 1 | treat == 2 | treat == 3)

###FIGURE 2#########

##helper function 

## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This does the summary. For each group's data frame, return a vector with
    # N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )

    # Rename the "mean" column    
    datac <- rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}


harass <- harass[!is.na(harass$SubLegStand),] #remove NA values from DV

data_sum1 <- summarySE(harass, measurevar="SubLegStand", groupvars=c("treat")) #get averages and CIs for substantive legitimacy across treatments 

# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.1) # move them .05 to the left and right

#plot values 

plot1<-ggplot(data_sum1, aes(x=1, y= SubLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= SubLegStand-se, ymax= SubLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=4, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Substantive Legitimacy \n (standardized within countries)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Women's rights issue (sexual harassment) \n Substantive legitimacy beliefs")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=12))+ 
    theme(plot.title = element_text(size = 12))


harass <- harass[!is.na(harass$ProLeg),]  #remove NA values from DV

data_sum2 <- summarySE(harass, measurevar="ProLegStand", groupvars=c("treat")) #get averages and CIs for prcedural legitimacy across treatments 


plot2<-ggplot(data_sum2, aes(x=1, y= ProLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= ProLegStand-se, ymax= ProLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=4, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Procedural Legitimacy \n (standardized within countries)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Women's rights issue (sexual harassment) \n Procedural legitimacy beliefs")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=12))+ 
    theme(plot.title = element_text(size = 12))

require(gridExtra)
# Set the PDF output file
pdf("Figure2.pdf", width = 8.2, height = 5)  # Width and height are in inches

grid.arrange(plot1 + annotate("text", x=0.9675, y=-0.05, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.28, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.25, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), plot2 + annotate("text", x= 0.9675, y=-0.42, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.45, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.38, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), ncol=2, nrow=1)


# Close the graphics device
dev.off()


####Figure 4 #### USA v. UK  


#subset data for UK and US

USA<-subset(harass, Country == "USA")
UK<-subset(harass, Country == "UK") 

data_sum1 <- summarySE(USA, measurevar="SubLegStand", groupvars=c("treat")) #get averages and CIs for substantive legitimacy across treatments 

# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.1) # move them .05 to the left and right

plot1<-ggplot(data_sum1, aes(x=1, y= SubLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= SubLegStand-se, ymax= SubLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=3, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Substantive Legitimacy \n (standardized within countries)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Women's rights issue (sexual harassment) \n Substantive legitimacy beliefs \n US sample")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=12))+ 
    theme(plot.title = element_text(size = 12))


harass <- harass[!is.na(harass$ProLeg),] #remove NA values from DV

data_sum2 <- summarySE(USA, measurevar="ProLegStand", groupvars=c("treat")) #get averages and CIs for procedural legitimacy across treatments 

plot2<-ggplot(data_sum2, aes(x=1, y= ProLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= ProLegStand-se, ymax= ProLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=3, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Procedural Legitimacy \n (standardized within countries)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Women's rights issue (sexual harassment) \n Procedural legitimacy beliefs \n US sample")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=12))+ 
    theme(plot.title = element_text(size = 12))


##UK

data_sum1 <- summarySE(UK, measurevar="SubLegStand", groupvars=c("treat")) #get averages and CIs for substantive legitimacy across treatments 


plot3<-ggplot(data_sum1, aes(x=1, y= SubLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= SubLegStand-se, ymax= SubLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=3, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Substantive Legitimacy \n (standardized within country)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Women's rights issue (sexual harassment) \n Substantive legitimacy beliefs \n UK sample")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=12))+ 
    theme(plot.title = element_text(size = 12))

data_sum2 <- summarySE(UK, measurevar="ProLegStand", groupvars=c("treat")) #get averages and CIs for procedural legitimacy across treatments 

plot4<-ggplot(data_sum2, aes(x=1, y= ProLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= ProLegStand-se, ymax= ProLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=3, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Procedural Legitimacy \n (standardized within country)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Women's rights issue (sexual harassment) \n Procedural legitimacy beliefs \n UK sample")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=12))+ 
    theme(plot.title = element_text(size = 12))


require(gridExtra)

# Set the PDF output file
pdf("Figure4.pdf", width = 8.6, height = 5.6)  # Width and height are in inches


grid.arrange(plot1 + annotate("text", x=0.9675, y=0.15, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.45, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.25, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), plot2 + annotate("text", x= 0.9675, y=-0.05, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.55, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.42, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), plot3 + annotate("text", x=0.9675, y=0.15, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.45, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.25, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), plot4 + annotate("text", x= 0.9675, y=-0.15, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.65, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.42, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), ncol=2, nrow=2)

# Close the graphics device
dev.off()


##FIGURE 5, non-women's rights issue 

##Subset treatments to respondents who saw animal mistreatment vignette

animal<-subset(MergedN, treat == 4 | treat == 5 | treat == 6)

animal <- animal[!is.na(animal $SubLegStand),] #remove NA values from DV

data_sum1 <- summarySE(animal, measurevar="SubLegStand", groupvars=c("treat")) #get averages and CIs for substantive legitimacy across treatments 


# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.1) # move them .05 to the left and right


plot3<-ggplot(data_sum1, aes(x=1, y= SubLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= SubLegStand-se, ymax= SubLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=4, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Substantive Legitimacy \n (standardized within countries)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Non women's rights issue (animal mistreatment) \n Substantive legitimacy beliefs")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=5))+ 
    theme(plot.title = element_text(size = 11))


animal <- animal[!is.na(animal $ProLeg),] #remove NA values from DV

data_sum2 <- summarySE(animal, measurevar="ProLegStand", groupvars=c("treat")) #get averages and CIs for procedural legitimacy across treatments 

plot4<-ggplot(data_sum2, aes(x=1, y= ProLegStand, group= treat)) + 
  geom_errorbar(aes(ymin= ProLegStand-se, ymax= ProLegStand +se), colour="black", width=.1, position=pd) +
    geom_line(position=pd) +
        geom_point(position=pd, size=4, shape=21, fill="white") + # 21 is filled circle
    xlab("") +
    ylab("Average Procedural Legitimacy \n (standardized within countries)") +
scale_y_continuous(limits = c(-1, 1)) +
    ggtitle("Non women's rights issue (animal mistreatment) \n Procedural legitimacy beliefs")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_text(size=5))+ 
    theme(plot.title = element_text(size = 11))

require(gridExtra)

pdf("Figure5.pdf", width = 8.6, height = 5.6)  # Width and height are in inches


grid.arrange(plot3 + annotate("text", x=0.9675, y=0.03, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.2, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.2, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), plot4 + annotate("text", x= 0.9675, y=-0.36, label= "All-Male Council", size=3.3) + annotate("text", x = 1, y=0.4, label = "Gender-Balanced \n Council", size=3.3) + annotate("text", x = 1.032, y=0.35, label = "Quota-Elected \n Gender-Balanced \n Council", size=3.3), ncol=2, nrow=1)

# Close the graphics device
dev.off()


##FOR SI

##TABLE SI.4 -- balance on women's issue treatments 

harass $Argentina <-recode(harass $Country, "'Argentina' = 1; else =0" ) 
harass $Australia <-recode(harass $Country, "'Australia' = 1; else =0" ) 
harass $Brazil <-recode(harass $Country, "'Brazil' = 1; else =0" ) 
harass $France <-recode(harass $Country, "'France' = 1; else =0" ) 
harass $Mexico <-recode(harass $Country, "'Mexico' = 1; else =0" ) 
harass $Norway <-recode(harass $Country, "'Norway' = 1; else =0" ) 
harass $NZ <-recode(harass $Country, "'NZ' = 1; else =0" ) 
harass $Peru <-recode(harass $Country, "'Peru' = 1; else =0" ) 
harass $Portugal <-recode(harass $Country, "'Portugal' = 1; else =0" ) 
harass $Spain <-recode(harass $Country, "'Spain' = 1; else =0" ) 
harass $UK <-recode(harass $Country, "'UK' = 1; else =0" ) 
harass $USA <-recode(harass $Country, "'USA' = 1; else =0" ) 


tgender<-tapply(harass $gender, harass $treat, mean, , na.rm=TRUE)
tideo<-tapply(harass $leftright_1, harass $treat, mean, , na.rm=TRUE)
targ<-tapply(harass $Argentina, harass $treat, mean, , na.rm=TRUE)
taus<-tapply(harass $Australia, harass $treat, mean, , na.rm=TRUE)
tbra<-tapply(harass $Brazil, harass $treat, mean, , na.rm=TRUE)
tfra<-tapply(harass $France, harass $treat, mean, , na.rm=TRUE)
tmex<-tapply(harass $Mexico, harass $treat, mean, , na.rm=TRUE)
tnor<-tapply(harass $Norway, harass $treat, mean, , na.rm=TRUE)
tnz<-tapply(harass $NZ, harass $treat, mean, , na.rm=TRUE)
tper<-tapply(harass $Peru, harass $treat, mean, , na.rm=TRUE)
tpor<-tapply(harass $Portugal, harass $treat, mean, , na.rm=TRUE)
tspa<-tapply(harass $Spain, harass $treat, mean, , na.rm=TRUE)
tuk<-tapply(harass $UK, harass $treat, mean, , na.rm=TRUE)
tusa<-tapply(harass $USA, harass $treat, mean, , na.rm=TRUE)


tapply_df <- as.data.frame(rbind(tgender, tideo, targ, taus, tbra, tfra, tmex, tnor, tnz, tper, tpor, tspa, tuk, tusa))

latex_table <-xtable(tapply_df)

# Export to a .tex file
sink("tableSI.4.tex")
print(latex_table, type = "latex")
sink()


##TABLE SI.5 -- balance on animal treatment treatments 

animal $Argentina <-recode(animal $Country, "'Argentina' = 1; else =0" ) 
animal $Australia <-recode(animal $Country, "'Australia' = 1; else =0" ) 
animal $Brazil <-recode(animal $Country, "'Brazil' = 1; else =0" ) 
animal $France <-recode(animal $Country, "'France' = 1; else =0" ) 
animal $Mexico <-recode(animal $Country, "'Mexico' = 1; else =0" ) 
animal $Norway <-recode(animal $Country, "'Norway' = 1; else =0" ) 
animal $NZ <-recode(animal $Country, "'NZ' = 1; else =0" ) 
animal $Peru <-recode(animal $Country, "'Peru' = 1; else =0" ) 
animal $Portugal <-recode(animal $Country, "'Portugal' = 1; else =0" ) 
animal $Spain <-recode(animal $Country, "'Spain' = 1; else =0" ) 
animal $UK <-recode(animal $Country, "'UK' = 1; else =0" ) 
animal $USA <-recode(animal $Country, "'USA' = 1; else =0" ) 


tgender<-tapply(animal $gender, animal $treat, mean, , na.rm=TRUE)
tideo<-tapply(animal $leftright_1, animal $treat, mean, , na.rm=TRUE)
targ<-tapply(animal $Argentina, animal $treat, mean, , na.rm=TRUE)
taus<-tapply(animal $Australia, animal $treat, mean, , na.rm=TRUE)
tbra<-tapply(animal $Brazil, animal $treat, mean, , na.rm=TRUE)
tfra<-tapply(animal $France, animal $treat, mean, , na.rm=TRUE)
tmex<-tapply(animal $Mexico, animal $treat, mean, , na.rm=TRUE)
tnor<-tapply(animal $Norway, animal $treat, mean, , na.rm=TRUE)
tnz<-tapply(animal $NZ, animal $treat, mean, , na.rm=TRUE)
tper<-tapply(animal $Peru, animal $treat, mean, , na.rm=TRUE)
tpor<-tapply(animal $Portugal, animal $treat, mean, , na.rm=TRUE)
tspa<-tapply(animal $Spain, animal $treat, mean, , na.rm=TRUE)
tuk<-tapply(animal $UK, animal $treat, mean, , na.rm=TRUE)
tusa<-tapply(animal $USA, animal $treat, mean, , na.rm=TRUE)


tapply_df <- as.data.frame(rbind(tgender, tideo, targ, taus, tbra, tfra, tmex, tnor, tnz, tper, tpor, tspa, tuk, tusa))

xtable(tapply_df)

latex_table <-xtable(tapply_df)

# Export to a .tex file
sink("tableSI.5.tex")
print(latex_table, type = "latex")
sink()



###SI Table SI.6: main cross-country treatment effects tables. Values can be found here: 

# Perform t-tests and extract relevant values
t1 <- t.test(SubLegStand ~ treat, data = harass[harass$treat %in% c(1, 2), ])
t2 <- t.test(SubLegStand ~ treat, data = harass[harass$treat %in% c(1, 3), ])
t3 <- t.test(ProLegStand ~ treat, data = harass[harass$treat %in% c(1, 2), ])
t4 <- t.test(ProLegStand ~ treat, data = harass[harass$treat %in% c(1, 3), ])

# Create a data frame with extracted values
results <- data.frame(
  Variable = c("SubLegStand", "SubLegStand", "ProLegStand", "ProLegStand"),
  Comparison = c("All Male vs. Gender-Balanced", 
                 "All Male vs. Quota Gender-Balanced", 
                 "All Male vs. Gender-Balanced", 
                 "All Male vs. Quota Gender-Balanced"),
  t_stat = c(t1$statistic, t2$statistic, t3$statistic, t4$statistic),
  p_value = c(t1$p.value, t2$p.value, t3$p.value, t4$p.value),
  CI_lower = c(t1$conf.int[1], t2$conf.int[1], t3$conf.int[1], t4$conf.int[1]),
  CI_upper = c(t1$conf.int[2], t2$conf.int[2], t3$conf.int[2], t4$conf.int[2])
)

# Assign row names based on test identifiers
rownames(results) <- c("T1: SubLegStand (1 vs. 2)", 
                       "T2: SubLegStand (1 vs. 3)", 
                       "T3: ProLegStand (1 vs. 2)", 
                       "T4: ProLegStand (1 vs. 3)")

# Convert the data frame to LaTeX table
latex_table <- xtable(results)

# Print LaTeX code
print(latex_table, type = "latex")

# Save to a .tex file
sink("tableSI.6.tex")
print(latex_table, type = "latex")
sink()


##SI Table SI.8: treatment effects in for the US and UK can be found here 


# Perform t-tests for the US
t1 <- t.test(SubLegStand ~ treat, data = USA[USA$treat %in% c(1, 2), ])
t2 <- t.test(SubLegStand ~ treat, data = USA[USA$treat %in% c(1, 3), ])
t3 <- t.test(ProLegStand ~ treat, data = USA[USA$treat %in% c(1, 2), ])
t4 <- t.test(ProLegStand ~ treat, data = USA[USA$treat %in% c(1, 3), ])

# Perform t-tests for the UK
t5 <- t.test(SubLegStand ~ treat, data = UK[UK$treat %in% c(1, 2), ])
t6 <- t.test(SubLegStand ~ treat, data = UK[UK$treat %in% c(1, 3), ])
t7 <- t.test(ProLegStand ~ treat, data = UK[UK$treat %in% c(1, 2), ])
t8 <- t.test(ProLegStand ~ treat, data = UK[UK$treat %in% c(1, 3), ])

# Create a data frame with extracted values
results <- data.frame(
  Country = c(rep("USA", 4), rep("UK", 4)),
  Variable = c(rep("SubLegStand", 2), rep("ProLegStand", 2),  # USA
               rep("SubLegStand", 2), rep("ProLegStand", 2)), # UK
  Comparison = c("All Male vs. Gender-Balanced", "All Male vs. Quota Gender-Balanced", 
                 "All Male vs. Gender-Balanced", "All Male vs. Quota Gender-Balanced", 
                 "All Male vs. Gender-Balanced", "All Male vs. Quota Gender-Balanced", 
                 "All Male vs. Gender-Balanced", "All Male vs. Quota Gender-Balanced"),
  t_stat = c(t1$statistic, t2$statistic, t3$statistic, t4$statistic, 
             t5$statistic, t6$statistic, t7$statistic, t8$statistic),
  p_value = c(t1$p.value, t2$p.value, t3$p.value, t4$p.value, 
              t5$p.value, t6$p.value, t7$p.value, t8$p.value),
  CI_lower = c(t1$conf.int[1], t2$conf.int[1], t3$conf.int[1], t4$conf.int[1], 
               t5$conf.int[1], t6$conf.int[1], t7$conf.int[1], t8$conf.int[1]),
  CI_upper = c(t1$conf.int[2], t2$conf.int[2], t3$conf.int[2], t4$conf.int[2], 
               t5$conf.int[2], t6$conf.int[2], t7$conf.int[2], t8$conf.int[2])
)

# Assign row names based on test identifiers
rownames(results) <- c("SI.8: USA - T1", "SI.8: USA - T2", "SI.8: USA - T3", "SI.8: USA - T4", 
                       "SI.8: UK - T5", "SI.8: UK - T6", "SI.8: UK - T7", "SI.8: UK - T8")

# Convert the data frame to LaTeX table
latex_table <- xtable(results,)

# Print LaTeX code
print(latex_table, type = "latex")

# Save to a .tex file
sink("tableSI.8.tex")
print(latex_table, type = "latex")
sink()


##SI Table SI.9: main treatment effects for animal mistreatment 

# Perform t-tests for animal mistreatment vignettes 
t1 <- t.test(SubLegStand ~ treat, data = animal[animal$treat %in% c(4, 5), ])
t2 <- t.test(SubLegStand ~ treat, data = animal[animal$treat %in% c(4, 6), ])
t3 <- t.test(ProLegStand ~ treat, data = animal[animal$treat %in% c(4, 5), ])
t4 <- t.test(ProLegStand ~ treat, data = animal[animal$treat %in% c(4, 6), ])

# Create a data frame with extracted values
results <- data.frame(
  Variable = c(rep("SubLegStand", 2), rep("ProLegStand", 2)),
  Comparison = c("All Male vs. Gender-Balanced", "All Male vs. Quota Gender-Balanced", 
                 "All Male vs. Gender-Balanced", "All Male vs. Quota Gender-Balanced"),
  t_stat = c(t1$statistic, t2$statistic, t3$statistic, t4$statistic),
  p_value = c(t1$p.value, t2$p.value, t3$p.value, t4$p.value),
  CI_lower = c(t1$conf.int[1], t2$conf.int[1], t3$conf.int[1], t4$conf.int[1]),
  CI_upper = c(t1$conf.int[2], t2$conf.int[2], t3$conf.int[2], t4$conf.int[2])
)

# Assign row names based on test identifiers
rownames(results) <- c("SI.9: Animal - T1", "SI.9: Animal - T2", "SI.9: Animal - T3", "SI.9: Animal - T4")

# Convert the data frame to LaTeX table
latex_table <- xtable(results)

# Print LaTeX code
print(latex_table, type = "latex")

# Save to a .tex file
sink("tableSI9.tex")
print(latex_table, type = "latex")
sink()


##For SI Table S.10 --

##model based specification with covariates: substantive legitimacy 

harassGB<-subset(harass, treat == 2 | treat == 3) #quota penatly
harassAM_GB<-subset(harass, treat == 1 | treat == 2)
harassAM_GBQ<-subset(harass, treat == 1 | treat == 3) 

harassAM_GBQ$treat<-recode(harassAM_GBQ$treat, "3=2") #so value just goes up by one, not two

#main models for outcome of substantive legitimacy 

model1<-lm(SubLegStand ~ treat+ gender + leftright_1 + Country, data= harassAM_GB) #b/n AMP and GB
model2<-lm(SubLegStand ~ treat+ gender + leftright_1 + Country,  data= harassAM_GBQ) #b/n AMP and GBQ
model3<-lm(SubLegStand ~ treat+ gender + leftright_1 + Country,  data= harassGB) #b/n GB and GBQ

screenreg(list(model1,  model2,  model3))
texreg(list(model1,  model2,  model3), digits=3)

# Run texreg and extract coefficients into a data frame
models <- list(model1, model2, model3)
texreg_output <- screenreg(models, digits = 3, include.ci = FALSE)

# Convert texreg output to a matrix format
extract_coefficients <- function(model) {
  summary_model <- summary(model)
  coefs <- summary_model$coefficients[, 1:2]  # Extract estimates & standard errors
  return(data.frame(Estimate = coefs[, 1], StdError = coefs[, 2]))
}

# Apply function to all models
model_results <- lapply(models, extract_coefficients)

# Combine results into a single table
coef_names <- rownames(model_results[[1]])
results_df <- data.frame(
  Coefficient = coef_names,
  Model1_Estimate = model_results[[1]]$Estimate,
  Model1_StdError = model_results[[1]]$StdError,
  Model2_Estimate = model_results[[2]]$Estimate,
  Model2_StdError = model_results[[2]]$StdError,
  Model3_Estimate = model_results[[3]]$Estimate,
  Model3_StdError = model_results[[3]]$StdError
)

# Convert data frame to xtable
latex_table <- xtable(results_df)

# Print LaTeX code
print(latex_table, type = "latex")

# Save to a .tex file
sink("tableSI10.tex")
print(latex_table, type = "latex")
sink()


##For SI Table S.11 --

##model based specification with covariates: procedural legitimacy

#main models for outcome of procedural legitimacy 

model1<-lm(ProLegStand ~ treat+ gender + leftright_1 + Country, data= harassAM_GB) #b/n AMP and GB
model2<-lm(ProLegStand ~ treat+ gender + leftright_1 + Country,  data= harassAM_GBQ) #b/n AMP and GBQ
model3<-lm(ProLegStand ~ treat+ gender + leftright_1 + Country,  data= harassGB) #b/n GB and GBQ

screenreg(list(model1,  model2,  model3))
texreg(list(model1,  model2,  model3), digits=3)

# Run texreg and extract coefficients into a data frame
models <- list(model1, model2, model3)
texreg_output <- screenreg(models, digits = 3, include.ci = FALSE)

# Convert texreg output to a matrix format
extract_coefficients <- function(model) {
  summary_model <- summary(model)
  coefs <- summary_model$coefficients[, 1:2]  # Extract estimates & standard errors
  return(data.frame(Estimate = coefs[, 1], StdError = coefs[, 2]))
}

# Apply function to all models
model_results <- lapply(models, extract_coefficients)

# Combine results into a single table
coef_names <- rownames(model_results[[1]])
results_df <- data.frame(
  Coefficient = coef_names,
  Model1_Estimate = model_results[[1]]$Estimate,
  Model1_StdError = model_results[[1]]$StdError,
  Model2_Estimate = model_results[[2]]$Estimate,
  Model2_StdError = model_results[[2]]$StdError,
  Model3_Estimate = model_results[[3]]$Estimate,
  Model3_StdError = model_results[[3]]$StdError
)

# Convert data frame to xtable
latex_table <- xtable(results_df)

# Print LaTeX code
print(latex_table, type = "latex")

# Save to a .tex file
sink("tableSI11.tex")
print(latex_table, type = "latex")
sink()

##For SI Table S.12: outcome substantive legitimacy 
##Robustness check, excluding those on the political left 

#excluding political liberals -- those that score five or above on the ideological self placement scale 

harassAM_GB_ModCon<-subset(harassAM_GB, leftright_1 > 4)
harassAM_GBQ_ModCon <-subset(harassAM_GBQ, leftright_1 > 4)
harassGB_ModCon <-subset(harassGB, leftright_1 > 4)

##models for substantive legitimacy

model1<-lm(SubLegStand ~ treat+ gender + leftright_1 + Country, data= harassAM_GB_ModCon) #b/n AMP and GB
model2<-lm(SubLegStand ~ treat+ gender + leftright_1 + Country,  data= harassAM_GBQ_ModCon) #b/n AMP and GBQ
model3<-lm(SubLegStand ~ treat+ gender + leftright_1 + Country,  data= harassGB_ModCon) #b/n GB and GBQ

screenreg(list(model1,  model2,  model3))
texreg(list(model1,  model2,  model3), digits=3)

# Run texreg and extract coefficients into a data frame
models <- list(model1, model2, model3)
texreg_output <- screenreg(models, digits = 3, include.ci = FALSE)

# Convert texreg output to a matrix format
extract_coefficients <- function(model) {
  summary_model <- summary(model)
  coefs <- summary_model$coefficients[, 1:2]  # Extract estimates & standard errors
  return(data.frame(Estimate = coefs[, 1], StdError = coefs[, 2]))
}

# Apply function to all models
model_results <- lapply(models, extract_coefficients)

# Combine results into a single table
coef_names <- rownames(model_results[[1]])
results_df <- data.frame(
  Coefficient = coef_names,
  Model1_Estimate = model_results[[1]]$Estimate,
  Model1_StdError = model_results[[1]]$StdError,
  Model2_Estimate = model_results[[2]]$Estimate,
  Model2_StdError = model_results[[2]]$StdError,
  Model3_Estimate = model_results[[3]]$Estimate,
  Model3_StdError = model_results[[3]]$StdError
)

# Convert data frame to xtable
latex_table <- xtable(results_df)

# Print LaTeX code
print(latex_table, type = "latex")

# Save to a .tex file
sink("tableSI12.tex")
print(latex_table, type = "latex")
sink()

##For SI Table S.13: outcome procerudal legitimacy 


model1<-lm(ProLegStand ~ treat+ gender + leftright_1 + Country, data= harassAM_GB_ModCon) #b/n AMP and GB
model2<-lm(ProLegStand ~ treat+ gender + leftright_1 + Country,  data= harassAM_GBQ_ModCon) #b/n AMP and GBQ
model3<-lm(ProLegStand ~ treat+ gender + leftright_1 + Country,  data= harassGB_ModCon) #b/n GB and GBQ

screenreg(list(model1,  model2,  model3))
texreg(list(model1,  model2,  model3), digits=3)

# Run texreg and extract coefficients into a data frame
models <- list(model1, model2, model3)
texreg_output <- screenreg(models, digits = 3, include.ci = FALSE)

# Convert texreg output to a matrix format
extract_coefficients <- function(model) {
  summary_model <- summary(model)
  coefs <- summary_model$coefficients[, 1:2]  # Extract estimates & standard errors
  return(data.frame(Estimate = coefs[, 1], StdError = coefs[, 2]))
}

# Apply function to all models
model_results <- lapply(models, extract_coefficients)

# Combine results into a single table
coef_names <- rownames(model_results[[1]])
results_df <- data.frame(
  Coefficient = coef_names,
  Model1_Estimate = model_results[[1]]$Estimate,
  Model1_StdError = model_results[[1]]$StdError,
  Model2_Estimate = model_results[[2]]$Estimate,
  Model2_StdError = model_results[[2]]$StdError,
  Model3_Estimate = model_results[[3]]$Estimate,
  Model3_StdError = model_results[[3]]$StdError
)

# Convert data frame to xtable
latex_table <- xtable(results_df)

# Print LaTeX code
print(latex_table, type = "latex")

# Save to a .tex file
sink("tableSI13.tex")
print(latex_table, type = "latex")
sink()

###SI Table SI.14 HTEs for substantive legitimacy 

#create a variable for those who identify as the most conservative: a self placement of 10

harassGB$Right<-recode(harassGB$leftright_1, "10 = 1; else =0" ) #sig for those 7, and

#add covar's
model1<-lm(SubLegStand ~ treat + gender + leftright_1 + Country, data= harassGB) #

#add interaction 1
model2<-lm(SubLegStand ~ treat + gender + leftright_1 + Country + I(gender*treat), data= harassGB) #quota penalty less for women

#add interaction 2
model3<-lm(SubLegStand ~ treat + gender + leftright_1 + Country + I(leftright_1 *treat), data= harassGB) #quota penalty less for women

model4<-lm(SubLegStand ~ treat + gender + Right + Country + I(treat* Right), data= harassGB) #quota penalty less for women

screenreg(list(model1,  model2,  model3, model4), stars = c(0.001, 0.01, 0.05, 0.1) )
texreg(list(model1,  model2,  model3, model4), digits=3, stars = c(0.001, 0.01, 0.05, 0.1) )



###SI Table SI.15 HTEs for procedural legitimacy 

#add covar's
model1<-lm(ProLegStand ~ treat + gender + leftright_1 + Country, data= harassGB) #

#add interaction 1
model2<-lm(ProLegStand ~ treat + gender + leftright_1 + Country + I(gender*treat), data= harassGB) #quota penalty less for women

#add interaction 2
model3<-lm(ProLegStand ~ treat + gender + leftright_1 + Country + I(leftright_1 *treat), data= harassGB) #quota penalty less for women

model4<-lm(ProLegStand ~ treat + gender + Right + Country + I(treat* Right), data= harassGB) #quota penalty less for women

screenreg(list(model1,  model2,  model3, model4), stars = c(0.001, 0.01, 0.05, 0.1) )
texreg(list(model1,  model2,  model3, model4), digits=3, stars = c(0.001, 0.01, 0.05, 0.1) )


